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Abstract. Floating point operations are fast, but require continuous effort on the part 
of the user in order to ensure that the results are correct. This burden can be shifted away 
from the user by providing a library of exact analysis in which the computer handles the 
error estimates. Previously, we |KS11] provided a fast implementation of the exact real 
numbers in the COQ proof assistant. Our implementation improved on an earlier imple- 
mentation by O'Connor |O'C08j by using type classes to describe an abstract specification 
of the underlying dense set from which the real numbers are built. In particular, we used 
dyadic rationals built from COQ's machine integers to obtain a 100 times speed up of the 
basic operations already. 

This article is a substantially expanded version of |KS11] in which the implementation 
is extended in the various ways. First, we implement and verify the sine and cosine 
function. Secondly, we create an additional implementation of the dense set based on COQ's 
fast rational numbers. Thirdly, we extend the hierarchy to capture order on undecidable 
structures, while it was limited to decidable structures before. This hierarchy, based on 
type classes, allows us to share theory on the naturals, integers, rationals, dyadics, and reals 
in a convenient way. Finally, we obtain another dramatic speed-up by avoiding evaluation 
of termination proofs at runtime. 



1. Introduction 

There is a big gap between numerical algorithms in research papers, which typically use 
concepts like Hilbert spaces and fixed point theorems from functional analysis, and their 
actual implementation, which uses floating point^ numbers and matrix operations. This gap 
makes it difficult to trust the code. Similarly, this gap is undesirable in proofs of theorems 
(e.g. Kepler conjecture [Hal02| . existence of the Lorentz attractor |Tuc02j ) that rely on 
the execution of this code for numerical approximation. Finally, from a purely software 
engineering point of view, the situation is undesirable, because the gap between the (abstract 

The research leading to these results has received funding from the European Union's 7th Framework 
Programme under grant agreement nr. 243847 (ForMath). 

^By floating points we mean numbers of the shape n*b'^, where n and e are integers with a finite precision 
and b is the base for scaling (typically 2, 10 or 16). The most widely used form of floating point arithmetic 
is the IEEE 754 standard, which is present in many hardware and software implementations. 
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mathematical) numerical algorithms and the (concrete floating point) implemented program 
makes the code difficult to maintain. 

The challenge to close this gap has already been posed by Bishop in his fundamental 
work on constructive analysis |Bis67) . Bishop proposed to use constructive analysis to 
bridge this gap. Moreover, we can narrow this gap by using 

• exact real numbers or intervals instead of floating point numbers; 

• functional programming instead of imperative programming; 

• dependent type theory which allows us to compute with complete metric spaces; 

• a proof assistant which allows us to verify the correctness proofs; 

• constructive mathematics to tightly connect mathematics with computations and 
to avoid computationally impossible case distinctions. 

Separately, all these tools have proved itself. By going to the limits of this proven technol- 
ogy we should be able to come within a small constant factor of floating point computa- 
tions. In this way one would obtain a tool suitable for research and education in numerical 
analysis that allows one to compute abstractly at the level of functional analysis, e.g. to 
compute fixed points of operators on Hilbert spaces. Like the development of Fortran 
and MATLAB this will require a huge amount of work. In the present paper we focus on 
the performance of real number computation in the COQ proof assistant. 

Real numbers, being infinite objects, cannot be represented exactly in a computer. 
Hence, in constructive analysis |Bis67j one uses functions which when fed a desired precision 
approximate a real numbers by a rational, or a dyadic number, to within that precision. 
The real numbers are the completion of the rationals. This completion construction can be 
organized in a monad, a familiar construct from functional programming. The completion 
monad provides an efficient combination of proving and computing [O'COTj . In this way, 
O'Connor |Q'C08j implements exact real numbers and the transcendental functions on them 
in COQ. A number of possible improvements in this implementation were already suggested 
in [QSTOllOTO] . 

(1) Use Coq's new machine integers instead of the integers built from ordinary inductive 
data types; 

(2) use dyadic rationals (that are numbers of the shape n * 2*^ for n,e G Z, also known 
as infinitary floats) instead of ordinary rationals; 

(3) use approximate division to improve the implementation of power series. 

Here we carry out all three optimizations. Unfortunately, changing O'Connor's implemen- 
tation to use the new machine integers was far from trivial, as he used a particular concrete 
representation of the rationals. To avoid this in the future, we provide an abstract spec- 
ification of the dense set as approximate rationals. Finally, we obtain another dramatic 
speed-up by avoiding evaluation of termination proofs at runtime. 

Outline. Section [2] describes some aspects of the COQ proof assistant relevant for our 
development. Section [3] describes metric spaces, monads, and O'Connor's implementa- 
tion of the real numbers |O'C07j . Section |4] extends Spitters and van der Weegen's ap- 
proach to abstract interfaces using type classes |SvdWl l] . Section [5] describes the the- 
ory of approximate rationals, our implementation of the real numbers, and deals with 
computing power series and the square root. We finish with some benchmarks in Sec- 
tion [6] and conclusions in Section [7j The sources of our developments can be found 
at [http : 7/ robbertkrebber s . nl/research/ reals| 
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2. The Coq-system 

The COQ proof assistant is based on tlie calculus of inductive constructions |CH88| 
ICP90j . a dependent type theory with (co)inductive types; see |Coq08[[BC04] . In true Curry- 
Howard fashion, it is both a pure functional programming language with an expressive type 
system, and a language for mathematical statements and proofs. We highlight some aspects 
of COQ relevant for our development. 

2.1. Types and propositions. Propositions in COQ are types |ML98|lML82j . which them- 
selves have types called sorts. COQ features a distinguished sort called Prop that one may 
choose to use as the sort for types representing propositions. The distinguishing feature of 
the Prop sort is that terms of non-Prop type may not depend on the values of inhabitants of 
Prop types (that is, proof terms). This regime of discrimination establishes a weak form of 
proof irrelevance, in that changing a proof can never affect the result of value computations. 
On a practical level, this lets COQ safely erase all Prop components when extracting certi- 
fied programs to OCAML or HASKELL. We should note however, that in practice, COQ's 
extraction mechanism |Let08| is still very hard to use for programs with the complexity, in 
terms of depth of definitions, that we are interested in |CFS03| ICFLOGj . 

2.2. Constructive indefinite description. Constructive indefinite description |BC04| 
14.2.3, 15.4] states that given a decidable predicate over the natural numbers, a Prop based 
existential can be converted into a Type based one. Its formal statement can be found in 
the standard library: 

Lemma constructive_indefinite_description_nat (P : nat -5> Prop) : 
(V X : nat, {P x} + P x}) ^ (3 n : nat, P n) ^ {n : nat | P n} 

Here the notation {x : A \ P x} for P : A — )• Prop denotes a S-type. This lemma can be 
seen as a form of Markov's principle in COQ. The algorithm does a bounded search for a 
new witness satisfying the predicate. The witness from the Prop based existential is only 
used to prove termination of the search. No information flows from the Prop universe to the 
Type universe because the witness found for the Type based existential is independent of the 
witness from the Prop based one. 

2.3. Equality, setoids, and rewriting. Because the COQ type theory lacks quotient types 
(as it endangers the decidability of type checking), one usually bases abstract structures 
on a setoid ('Bishop set'): a type equipped with an equivalence relation [Bis67| IHof97] . 
This leads to a naive set theory as described by Palmgren |Pal09| . When the user attempts 
to substitute a given (sub)term using an equality, the system keeps track of, resolves, and 
combines proofs of equivalence [Soz09] . 

The 'native' notion of equality in COQ, Leibniz equality, is that of terms being convert- 
ible, naturally reified as a proposition by the inductive type family eq with single constructor 
eq_refl : V (T : Type)(x : T), x = x, where a = b is notation for eq T a b. Since convertibility is a con- 
gruence, a proof of a = b lets us substitute b for a anywhere inside a term without further 
conditions. Our interest is in more complicated equalities, so we diverge from COQ tradition 
and reserve = for setoid equality. Rewriting with = does give rise to side conditions. For 
instance, consider formal fractions of integers as a representation of rationals. Rewriting a 
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subterm using such an equality is permitted only if the subterm is an argument of a func- 
tion that has been proven to respect the equality. Such a function is called Proper, and that 
property must be proved for each function in whose arguments we wish to enable rewriting. 

2.4. Type classes. Type classes have been a great success story in the Haskell functional 
programming language, as a means of organizing interfaces of abstract structures. COQ's 
type classes provide a superset of their functionality, but are implemented in a different 
way. 

In Haskell and Isabelle, type classes and their instances are second class. They 
are handled as specialized syntactic constructs whose semantics are given specifically by 
the type class apparatus. By contrast, the expressivity of dependent types and inductive 
families as supported in COQ, combined with the use of pre-existing technology in the 
system (namely proof search and implicit arguments) enable a first class type class imple- 
mentation [SO08J: classes are ordinary record types ('dictionaries'), instances are ordinary 
constants of these record types (registered as hints with the proof search machinery), class 
constraints are ordinary implicit parameters, and instance resolution is achieved by aug- 
menting the unification algorithm to invoke ordinary proof search for implicit arguments of 
class type. Thus, type classes in COQ are realized by relatively minor syntactic aids that 
bring together existing facilities of the theory and the system into a coherent idiom, rather 
than by introduction of a new category of qualitatively different definitions with their own 
dedicated semantics. 

We use the algebraic hierarchy based on type classes and its abstract specification of 
N, Z and Q described in [SvdWll] . Unfortunately, we should note that we have clearly met 
the efficiency problems connected to the current implementation of type classes in COQ. 
Luckily, these efficiency problems are limited to instance resolution which is only performed 
at compile time. Type classes have only a very minor eff'ect on the computation time of 
type checked terms due to the absence of code inlining; see Section [6] for timings. 

2.5. Virtual machine and machine integers. COQ includes a virtual machine |GL02j . 
vm.compute, based on Ocaml's virtual machine to allow efficient evaluation. Both the ab- 
stract machine and its compilation scheme have been proved correct, in COQ, with respect 
to the weak reduction semantics. However, we still need to extend our trusted core to a 
bigger kernel, as the implementation has not been formally verified. 

Machine integers were also added to the COQ system [AGSTIO] . The usual evaluation 
inside COQ (compute) uses a special inductive type for cyclic integers, but the virtual machine 
uses Ocaml's machine integers. This allows for a big speed-up, for which we pay by having 
to trust (the virtual machine and) that Ocaml treats these integers correctly. The time 
diff'erence between computation with COQ's machine integers and Ocaml's BigJnt is about 
a factor of 20 |Spill| on primality tests. 

3. Metric spaces 

Having completed our brief description of the COQ-system, we now turn to O'Connor's 
formalization of exact real numbers |O'C07j . Traditionally, a metric space is defined as a set 
X with a metric function d : X x X ^ IR"*" satisfying certain axioms. We use a more relaxed 
definition of a metric space that does not require the metric be a function; see also [RicOSj . 
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The metric is represented via a (respectful) ball relation B : — )■ X — )■ X — )• Prop 
satisfying: 

msp.refl : Vxe, a; a; 

msp_sym :Wxye, Ji^x y y x 

mspJriangle '.Wxy z ei e2, B^^ x y ^ B^^ y z ^ B^^^^^ x z 
msp_closed :\lxye, (V5, B^^^ x y) ^ 'B^x y 
msp_eq : "ixy, (Ve, Ji^x y) x = y 

The ball relation x y expresses that the points x and y are within e of each other. We 
call this a ball relationship because the partially applied relation B^ x : X ^ Prop is a 
predicate that represents the closed ball of radius e around the point x. For example, the 
ball relation on Q is B^ x y := \x — y\ < e. 
A metric space X is a prelength space if: 

\/abe6i62, e < 61 + 62 ^ B^a b ^ 3c, B^^ a c A Bg^cb. 

This property states that if two points a and b within e of each other, then there exists 
curve of length d(a, b) in the completion of X that connects a and b. The metric space Q 
is a prelength space. 

We will introduce the completion of a metric space as a monad. In order to do this we 
will first introduce monads. 



3.1. Monads. Moggi Mog89 recognized that many non-standard forms of computation 
may be modeled by monads'^. Wadler |Wad92] popularized their use in functional program- 
ming. Monads are now an established tool to structure computation with side-effects. For 
instance, programs with input X and output Y which have access to a mutable state S can 
be modeled as functions of type XxS^YxS, or equivalently X ^ {Y x S)''^ . The type 
constructor VJIY := (Y x S)''^ is an example of a monad. Similarly, partial functions may 
be modeled by maps X ^ Y±, where Yl := y -|- () is a monad. 

The formal definition of a (strong) monad is a triple (9Jt, return, bind) consisting of a 
type constructor Tl and two functions: 

return : X MX 

bind : {X mY) TlX MY 

We will denote return x as x, and bind / as /. These two operations must satisfy: 

bind return a = a 
fa = fa 
f{ga) = bind (fog) a 



in category theory one would speak about the Kleish category of a (strong) monad. 
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3.2. Completion monad. The completion of a metric space X is defined by: 

:= {/ : Q+ ^ X I Vei^s, B,^+,^ (/ei) {fe^)]. 

Given metric spaces X and Y , a function f : X \s uniformly continuous with modulus 
fj,f:Q+^ Q+ if: 

VexiX2, B^^.^xi X2 ^B^{f xi) (/X2). 
Completion is a monad on the category of metric spaces with uniformly continuous func- 
tions. The function return : X — )• <tX defined by Xxe, x is the embedding of a metric space 
in its completion. Moreover, a uniformly continuous function / : X — )• ^Y with modulus 
fj,f can be lifted to operate on complete metric spaces as bind / : (tX — )■ (tY defined by 
Ax e, / (x (;U/|)) |. A restriction to prelength spaces is essential for this efficient definition 
of bind. 

One advantage of this approach is that it helps us to work with simple representations. 
Let IR := CQ. Then to specify a function from IR, — )• R, we define a uniformly continuous 
function / : Q — )• R, and obtain / : IR, — )• IR as the required function. Hence, the comple- 
tion monad allows us to do in a structured way what was already folklore in constructive 
mathematics: to work with simple, often decidable, approximations to continuous objects. 

4. Abstract interfaces using type classes 

An important part of this work is to further develop the algebraic hierarchy based on 
type classes by Spitters and van der Weegen |SvdWl l] . Especially, we extend their hierarchy 
with constructive fields, order theory and interfaces for mathematical operations, such as 
shift and power, common in programming languages. This layer of abstraction makes both 
proof engineering and programming more flexible: it avoids duplication of code, it introduces 
a canonical way to refer to operations and properties, both by names and notations, and 
it allows us to easily swap different implementations of number representations and their 
operations. First we will briefly recap the design decisions made in [SvdWll] . 

Algebraic structures are expressed in terms of a number of carrier sets, a number of 
relations and operations, and a number of laws that the operations satisfy. One way of 
describing such a structure is by a bundled representation: one uses a dependently typed 
record that contains the carrier, operations and laws. For example a semigroup can be 
represented as follows. (The fields sg_car and sg_proper support our explicit handling of naive 
set theory in type theory.) 

Record SemiGroup : Type := { 
sg_car :> Setoid ; 
sg_op : sg_car sg_car sg_car ; 
sg_proper : Proper ((=) ^ (=) ^ (=)) sg_op ; 
sg_ass : V X y z, sg_op x (sg_op y z) = sg_op (sg_op x y) z) } 

However, this approach has some serious limitations, the most important one being a 
lack of support for sharing components. For example, suppose we group together two 
CommutativeMonoids in order to create a SemiRing. Now awkward hacks are necessary to es- 
tablish equality between the carriers. A second problem is that if we stack up these records 
to represent higher structures the projection paths become increasingly long. 

Historically these problems have been an acceptable trade-off because unbundled rep- 
resentations, in which the carrier and operations are parameterized, introduce even more 
problems. 
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Record SemiGroup {A} (e : A ^ A Prop) (sg_op : A ^ A ^ A) : Prop := { 
sg_proper : Proper (e e => e) sg_op ; 
sg_ass : V X y z, e (sg_op x (sg_op y z)) (sg_op (sg_op x y) z) } 

There is nothing to bind notation to, no structure inference, and declaring and passing 
requires too much manual bookkeeping. Spitters and van der Weegen have proposed a 
use of Coq's new type class machinery that resolves many of the problems of unbundled 
representations. Our current experiment confirms that this is a viable approach. 

An alternative solution is provided by packed classes [GGMROO] which use an alter- 
native, and older, implementation of a semblance of type classes: canonical structures; see 
also Section [7j Yet another approach would be to use modules. However, as these are not 
fist class, we would be unable to define, e.g. homomorphisms between algebraic structures. 

An operational type class is defined for each operation and relation. 

Class Equiv A := equiv: relation A. 
Infix "=" := equiv: type_scope. 
Class RingPlus A := ring_plus: A — > A — > A. 
Infix "+" := ring.plus. 

Now an algebraic structure is just a type class living in Prop that is parametrized by its 
carrier, relations and operations. This class contains all laws that the operations should 
satisfy. Since the operations are unbundled we can easily support sharing. For example let 
us consider the SemiRing interface. 

Class SemiRing A {e : Equiv A} {plus: RingPlus A} 

{mult: RinglVlult A} {zero: RingZero A} {one: RingOne A} : Prop := { 
semiring_mult_monoid :> @CommutativeMonoid A e mult one ; 
semiring_plus_monoid :> ©CommutativelVlonoid A e plus zero ; 
semiring.distr :> Distribute (.*.) (+) ; 
semiring_left_absorb :> LeftAbsorb (.*.) }. 

Without type classes it would be a burden to manually carry around the carrier, relations 
and operations. However, because these parameters are just type class instances, the type 
class machinery will perform that job for us. For example, 
Lemma example '{SemiRing R} x : 1 * x = x + 0. 

The backtick instructs COQ to automatically insert implicit declarations, namely e plus 
mult zero one. It also lets us omit a name for the SemiRing R parameter itself. All of these 
parameters will be given automatically generated names that we will never refer to. Further- 
more, instance resolution will automatically find instances of the operational type classes 
for the written notations. Thus the above is really: 

Lemma example {R e plus mult zero one} {P : ©SemiRing R e plus mult zero one} x : 
@equiv R e 
(@ring_mult R mult (@ring_one R one) x) 
(@ring_plus R plus x (@ring_zero R zero)). 

The syntax :> in the definition of SemiRing declares certain fields as substructures. This 
syntax should not be confused with the similar syntax for coercions in records (e.g. in the 
bundled representation of a SemiGroup on pagejo]). In current versions of COQ, inference 
of substructures is based on backward reasoning. In this example that means, each time a 
CommutativeMonoid R instance is needed, instance search may try to find a SemiRing R instance. 
This style of instance search presents some problems, as the following example illustates. 
Class Setoid_IVIorphism {A B} {Ae: Equiv A} {Be: Equiv B} (f: A ^ B) := { 
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setoidmor_a :> Setoid A; 
setoidmor_b :> Setoid B; 
sm_proper :> Proper ((=) (=)) f }. 

Each time we have to estabhsh Setoid R, instance search might try to find a Setoid.lVlorphism 
to or from R. Since this quickly results in a serious blow-up, we omit the :> declaration. 
Support for forward reasoning would solve this problem. If we would be in a context in 
which we know something to be a Setoid.lViorphism, then forward reasoning automatically 
infers that the source and target are Setoids. Ongoing work by Matthieu Sozeau should 
make this style of instance search available in COQ. 

Proving that an actual structure is an instance of the SemiRIng interface is straightfor- 
ward. First we define instances of the operational type classes. 

Instance nat.equiv: Equiv nat := eq. 
Instance nat_plus: RingPlus nat := plus. 
Instance nat_0: RingZero nat := 0%nat. 
Instance nat_1 : RingOne nat := 1%nat. 
Instance nat_mult: RingMult nat := mult. 

Here we see that instances are just ordinary constants of the class types. However, we use 
the Instance keyword instead of Definition to let the type class machinery register the instance. 
Now, to prove that the Peano naturals are in fact a semiring, we just write: 

Instance: SemiRIng nat. 
Proof. . . . Qed. 

This approach to interfaces proved useful to formalize a standard algebraic hierarchy. 
Combined with category theory and universal algebra, N and Z are represented as interfaces 
specifying an initial semiring and initial ring jSvdWll] . 

Class NaturalsToSemlRIng (A : Type) := 

naturals_to_semiring: V B '{RingMult B} '{RingPlus B} '{RingOne B} '{RingZero B}, A -s> B. 
Class Naturals A {e plus mult zero one} '{U: NaturalsToSemlRIng A} := { 

naturals.ring :> @SemiRing A e plus mult zero one; 

naturals_to_semiring_mor :> V '{SemlRing B}, SemiRing_Morphism (naturals_to_semiring A B); 
naturalsJnitial :> Initial (semirings. object A) }. 

These abstract interfaces for the naturals and integers make it easier to change the concrete 
representation in the future. No such simple specification for Q seems to exists, so we 
choose to specify it as the field of fractions of Z. More precisely, Q is specified as a field 
containing Z that moreover can be embedded into the field of fractions of Z. 

Inductive Frac R '{e : Equiv R} '{zero : RingZero R} : Type := 

frac { num : R ; den : R ; den.nonzero : den / }. 
Class RationalsToFrac (A : Type) := rationals_to_frac : V B '{Integers B}, A -s> Frac B. 
Class Rationals A {e plus mult zero one opp inv} '{U : ! RationalsToFrac A} : Prop := { 

rationalsJield :> @DecField A e plus mult zero one opp inv ; 

rationalsJrac :> V '{Integers Z}, Injective (rationals_to_frac A Z) ; 

rationals_frac_mor :> V '{Integers Z}, SemiRlng.lVlorphlsm (rationalsJoJrac A Z) ; 

rationals.embedJnts :> V '{Integers Z}, Injective {integers_to_ring Z A) }. 

4.1. Constructive fields and apartness. In constructive mathematics, the common no- 
tion of inequality as the negation of equality is often too weak because a proof of a negation 
lacks computational content. For example, in order to define the reciprocal of x G R, one 
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needs a witness e G Q+ that |x| > e. Such a witness cannot be extracted from a proof of 
X / 0. To solve this problem, one uses a setoid equipped with an apartness (irreflexive, 
asymmetric and co-transitive) relation describing inequality |TvD88j . 

The algebraic hierarchy in the CoRN library |CFGW04] has been built on top of such 
setoids. Unfortunately, this hierarchy is quite 'heavy' in practice. First, for structures with 
decidable equality, the negation of equality is the only tight apartness. Hence, when working 
with decidable structures, an apartness relation is unnecessary. Secondly, CoRN uses an 
informative (that is. Type based) apartness relation to facilitate extraction of witnesses. 
However, Coq's present implementation of setoid rewriting does not support rewriting over 
relations in Type. So, it does not allow us to replace equations in expressions involving 
CoRN's informative apartness and thus many proofs involve a lot of manual labor. 

To remedy these issues we propose an alternative solution. We use a non-informative 
(that is, Prop-based) apartness relation to enable setoid rewriting and include it just in the 
parts of the algebraic hierarchy where we actually need it. The latter keeps our interfaces 
clean and easy to use and should combine the best of two worlds. Type classes are of great 
help to reduce bookkeeping and clutter in proofs. 

Although using a non-informative apartness relation enables setoid rewriting, it dis- 
ables extraction of witnesses. Fortunately, in case of the reals, a witness can be obtained 



inefficiently by bounded linear search (see Section 2.2 and 5.1). We think our approach is 



a reasonable trade-off since the amount of reasoning exceeds the potential use of apartness 
in computation. In case we need a witness for efficient computation, we just have to spec- 
ify it explicitly. This approach of specifying witnesses explicitly was already preferred by 
O'Connor |O'C08j . even when an informative apartness was available. 

Our interface for a setoid with apartness (henceforth StrongSetoid) is as follows. 
Class Apart A := apart: relation A. 

Infix "X " := apart (at level 70, no associativity) : type.scope. 



Class StrongSetoid A {e: Equiv A} '{ap : Apart A} : Prop := { 
strong_setoid_irreflexive :> Irreflexive (x) ; 
strong_setoid_symmetric :> Symmetric (x) ; 
strong_setoid_cotrans :> CoTransitive (x) ; 
tight_apart :Vxy, ^xxyox = y}. 

This interface is equipped with a tight equality. We prove that each StrongSetoid is a Setoid. 
For decidable structures, we define the following class to describe that the apartness relation 
is the negation of equality. 

Class TrivialApart R '{Equiv R} {ap : Apart R} := trivial_apart :Vxy, xxyox/y. 

Given a setoid with decidable equality we can easily extend it to a StrongSetoid. 

Instance default-apart '{Equiv A} : Apart A | 20 := (/ ). 
Instance default_apart_trivial '{Equiv A} : TrivialApart A (ap:=default_apart). 
Lemma dec_strong_setoid '{Setoid A} '{Apart A} 
'{ITrivialApart A} '{V x y, Decision (x = y)} : StrongSetoid A. 

Unfortunately, the type class mechanism is unable to detect simple loops. Hence we define 
the above as an ordinary Lemma instead of an Instance. This trick prevents COQ from using it 
in instance search and therefore avoids endless derivations of the form StrongSetoid A, Setoid A, 
StrongSetoid A, . . . 
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For ordinary setoids we want functions to be Proper, which means that they respect 
equahty. For setoids with apartness we need a stronger property, strong extensionality. 

Class StrongSetoid.Morphism {A B : Type} {Ae: Equiv A} {Aap: Apart A} 
{Be: Equiv B} {Bap : Apart B} (f: A ^ B):= { 
strong_setoidmor_a: StrongSetoid A ; 
strong_setoidmor_b: StrongSetoid B ; 
strong.extensionality :Vxy, fxxfy-^xxy}. 

We prove that for each StrongSetoid_Morphism f we have Proper ((=)=> (=))f. The only structures 
for which we actuahy need apartness are implementations of the real numbers, hence we 
only base the Field class on top of a StrongSetoid instead of the complete algebraic hierarchy. 
Our class for fields is as follows. (The PropHolds class is explained in the next subsection.) 

Class Field A {e plus mult zero one opp} {ap: Apart A} {multJnv: Multlnv A} : Prop := { 
field_ring:> @Ring A e plus mult zero one opp ; 
field_strongsetoid:> StrongSetoid A ; 
field_plus_ext:> StrongSetoid.BinarylVlorphism (+) ; 
field_mult_ext:> StrongSetoid.BinarylVlorphism (.*.) ; 
field_nontrivial:> PropHolds (1 >c 0) ; 
mult_inv_proper:> Setoid_Morphism (//) ; 
mult_inverse: V x, 'x // x = 1 }. 

We do not include strong extensionality of the inverse and the reciprocal because it can 
be derived. For convenience, we define a class for fields with decidable equality whose 
reciprocal function is total. This class integrates nicely with COQ's rational numbers Q and 
bigQ, and the field tactic. This total reciprocal function should satisfy /O = 0, so properties 
as f{/x) = /{fx), /(/x) = X and /x * /y = /{x * y) hold without any additional premises. 
A diagram of our complete algebraic hierarchy is displayed in Figure [T] 

4.2. Order theory. Previous COQ libraries for ordered algebraic structures turn out to 
be too limited to abstract from N, Z, Q and IR and their various implementations. The 
formalization of ordered fields in the CoRN library |CFGW04] restricts to a very specific part 
of the algebraic hierarchy (namely fields). Letouzey's Numbers library, which is included 
in recent versions of COQ trunk, only considers N and Z. The Ssreflect library presently 
restricts to decidable structures with Leibniz equality. Moreover, even mathematically, the 
most convenient abstraction is not entirely clear. Lesnik [LeslOj provides a smooth order 
theoretic characterization of these structures as so-called streaks. We, however, prefer our 
theory below as it avoids partial functions. 

In this work we present a library that captures the notion of order on a variety of 
structures, including structures with undecidable equality. One of the building blocks of 
our hierarchy is a pseudo order, which is the constructive variant of a total order. 

Class PseudoOrder '{e : Equiv A} '{ap : Apart A} (so : Lt A) : Prop := { 
pseudo_order_setoid : StrongSetoid A ; 
pseudo_order_antisym : V x y, ^ (x < y a y < x) ; 
pseudo_order_cotrans :> CoTransitive {<) ; 
apart Jff_total_lt :Vxy, xxyox<yvy<x}. 

In case equality is decidable, this interface is rather awkward to work with. Therefore we 
present ways to go back and forth between the usual classical notions and their constructive 
variants. For example, we use the type class machinery to infer the classical trichotomy 
property. 



TYPE CLASSES FOR EFFICIENT EXACT REAL ARITHMETIC IN COQ 



11 




(a) The algebraic hierarchy (b) The order hierarchy 



Figure 1. The algebraic and order hierarchy. Dotted hnes denote derived 
inheritance, fihed nodes denote presence of apartness. 

Instance ItJrichotomy '{PseudoOrder A} '{ITrivialApart A} '{V x y, Decision (x = y)} : Trichotomy (<). 

Also, we can go the other way around. If we have a StrictSetoidOrder (an ordinary strict order 
built upon a setoid) satisfying the trichotomy property, we obtain a pseudo order. 

Lemma dec_strict_pseudo_order '{StrictSetoidOrder A} '{Apart A} 

'{ITrivialApart A} '{V x y, Decision (x = y)} '{ITrichotomy (<)}: PseudoOrder {<). 

In order to avoid loops, we define the above as an ordinary Lemma instead of an Instance. 
Next, one could extend a pseudo order to the standard notion of a (pseudo) ring order. 

Class PseudoRingOrder '{Equiv A} '{Apart A} '{RingPlus A} '{RingMult A} '{RingZero A} (Alt : Lt A) := { 
pseudo_ringorder_spo :> PseudoOrder Alt ; 
pseudo_ringorder_mult_ext :> StrongSetoid.BinaryMorphism (.*.) ; 
pseudo_ringorder_plus :> V z, StrictlyOrderPreserving (z +) ; 
pseudo_ringorder_mult :Vxy, 0<x-^0<y-^0<x*y}. 

However, we would like to apply our order library to implementations of the naturals, which 
are merely semirings, too. Therefore, we strengthen the above with a partial subtraction 
function (living in Prop, because we never use it for computations) and require addition to 
be order reflecting. We call this, apparently new notion, a PseudoSemiRingOrder. 

Class PseudoSemiRingOrder '{Equiv A} '{Apart A} '{RingPlus A} 
'{RingMult A} '{RingZero A} (Alt : Lt A) := { 
pseudo_srorder_strict :> PseudoOrder Alt ; 
pseudo_srorder_partiaLminus ;Vxy, ^y<x^3z, y = x + z; 
pseudo_srorder_plus :> V z, StrictOrderEmbedding (z +) ; 
pseudo_srorder_mult_ext :> StrongSetoid_BinarylVlorphism (.*.) ; 
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pseuclo_srorder_pos_mult_compat :Vxy, 0<x->0<y^O<x*y}. 

Instead of including the PseudoRingOrder class in our development, we include a lemma to 
construct a PseudoSemiRingOrder from a ring satisfying the PseudoRingOrder axioms. 

Given a pseudo (semiring) order, one could define the non-strict order x < y in terms of 
the strict order, namely as ^ y < x. However, this is quite inconvenient in practice, because 
we also want to talk about a priori different non-strict orders such as those defined in the 
standard library. Hence we introduce the following class. 

Class FullPseudoSemiRingOrder '{Equiv A} '{Apart A} '{RingPlus A} 
'{RingMult A} '{RingZero A} (Ale : Le A) (Alt ; Lt A) := { 
full_pseudo_srorder_pso :> PseudoSemiRingOrder Alt ; 
full_pseudo_srorderJe_iff_notJt_flip :Vxy, x<yo^y<x}. 

A diagram of our complete order hierarchy is displayed in Figure [l] 

Our theory on abstract orders avoids duplication of theorems and proofs. For example, 
the following lemmas apply to N, Z, Q and the dyadics, because all of these structures 
form a FullPseudoSemiRingOrder. 

Lemma plus_compat Xi yi X2 y2 : Xi < yi X2 < y2 ^ Xi + X2 < yi + y2 • 

Lemma lt_1 _2 : 1 < 2. 

Lemma square.nonneg x : < x * x. 

To allow us to refer by canonical names to common properties, we introduce classes like 
those shown below: 

Class OrderPreserving {A B} {Ae: Equiv A} {Ale: Le A} {Be: Equiv B} {Ble: Le B} (f : A ^ B) := { 

order_preserving_morphism :> Order.Morphism ; 

order_preserving : '(x < y — ^ f x < f y) }■ 
Class OrderReflecting {A B} {Ae: Equiv A} {Ale: Le A} {Be: Equiv B} {Ble: Le B} (f : A ^ B) := { 

order_preserving_back_morphism :> Order.Morphism ; 

order_preserving_back : '(f x < f y x < y) }. 

Here, an Order_Morphism is just the factoring out of the common parts of both classes; namely 
that f and < respect equality. For the case of multiplication these properties have additional 
premises, for example: 

Global Instance: V (z : R), PropHolds (0 < z) OrderPreserving (z *.). 

We introduce the PropHolds class to let the type class machinery prove these properties 
automatically. For example consider: 

Lemma example (n : N) (x y : R) : x < y — > (2 " n + 2) * x < (2 " n + 2) * y. 
Proof, intros. now apply (order_preserving (2 " n + 2 .*)). Qed. 

In order to use order_preserving, we need a proof of PropHolds (0 < 2 " n + 2). Type class resolu- 
tion is able to prove this in a fully automated way because we have the following instances: 
Instance: PropHolds (0 < 2); 

Instance: V x y : R, PropHolds (0 < x) ^ PropHolds (0 < y) ^ PropHolds (0 < x + y) 
Instance: V (n : N) (x : R), PropHolds (0 < x) ^ PropHolds (0 < x " n) 

For arbitrary instances of N, Z, Q it is easy to define an order satisfying these interfaces: 

Instance natJe '{Naturals N} : Le N | 10 := A x y, 3 z, y = x + z. 
Instance natJt '{Naturals N} : Lt N | 10 := A x y, x < y A x / y. 

However, often we encounter an a priori different order on a structure, most likely an 
order defined in Coq's standard library (like NIe and Nit on N). Therefore we prove that a 
FullPseudoSemiRingOrder uniquely specifies the order on N, Z and Q. For example: 
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Context '{Naturals N} '{Naturals N2} {f : N ^ N2} '{!SemiRing_Morphism f} 
'{Apart N} '{ITrivialApart N} '{IFullPseudoSemlRingOrder (A:=N) NIe Nit} 
'{Apart N2} '{ITrivialApart N2} '{IFullPseudoSemlRingOrder (A;=N2) N2le N2lt}. 

Global Instance: OrderEmbedding f. 

Unfortunately COQ has no support to have an argument be 'inferred if possible, generalized 
otherwise'; see [SvdWll] . When declaring a parameter of FuilPseudoSemiRlngOrder, one is 
often in a context where most of its components are already available. Usually, only the 
parameters Le, Lt and Apart have to be introduced. The current workaround in these cases 
involves providing names for components that are then never referred to, which is a bit 
awkward. In the above it would much nicer to write: 

Context '{Naturals N} '{Naturals N2} {f : N ^ N2} '{!SemiRing_Morphism f} 

'{ITrivialApart N} '{IFullPseudoSemiRlngOrder N} '{ITrivialApart N2} '{IFullPseudoSemlRingOrder N2}. 
Global Instance: OrderEmbedding f. 

4.3. Basic operations. The operation nat_pow is most commonly, but inefficiently, defined 
as repeated multiplication and the operation shift! is defined as repeated duplication. Instead 
we specify the desired behavior of these operations. This approach allows for different imple- 
mentations for different number representations and avoids definitions and proofs becoming 
implementation dependent. 

We introduce interfaces that specify the behavior of the operations abs, shifti, nat_pow and 
int_pow. Again there are various ways of specifying these interfaces: with S-types, bundled 
or unbundled. In general, S-types are convenient for functions whose specification is easy, 
for example: 

Class Abs A '{Equiv A} '{Le A} '{RingZero A} '{Grouplnv A} 

:= abs_sig: V (x : A), { y : A | (0 < x ^ y = x) a (x < ^ y = -x)}. 
Definition abs '{Abs A} := A x : A, ' (abs.sig x). 

However, for more complex operations, such as shifti, such an interface is different from the 
usual mathematical specification because we cannot quantify over all possible input values. 
Now there are two ways: a bundled or an unbundled interface. Since these interfaces are 
not used for hierarchies the disadvantages of the former do not apply. Let us first describe 
the former approach. 

Class ShiftL A B '{Equiv A} '{Equiv B} '{RingOne A} 

'{RingPlus A} '{RingMult A} '{RingZero B} '{RingOne B} '{RingPlus B} := { 
shifti : A ^ B ^ A ; 

shiftLproper : Proper {(=) =^ (=) ^ (=)) shifti ; 
shiftLO :> Rightldentity shifti ; 
shiftLS : V X n, shifti x (1 + n) = 2 * shifti x n }. 
Infix "< " := shifti (at level 33, left associativity). 

Although this interface seems reasonable, it does not work well in COQ. The simpi tactic 
which is used to simplify a goal will unfold occurrences of shifti to their underlying definition 
(for example in case of BigN, the expression x<n becomes BigN. shifti x n). This is rather 
inconvenient because COQ will then be unable to use lemmas concerning ^ for rewriting. 
This problem is caused because shifti is a projection of a record, which is in fact an t-redex 
(reduction of pattern-matching over a constructed term) that will be unfolded by simpl. Cur- 
rently there seems to be no way to adjust the behavior of simpl to remove this inconvenience. 
A similar problem was already observed in Ssreflect [GMTOS] . 
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Instead we use an unbundled interface, which has a lot in common with our interfaces 
for algebraic structures. Now shifti no longer contains an i-redex. 
Class ShiftL A B := shifti: A -> B -> A. 
Infix "< " := shifti (at level 33, left associativity). 

Class ShiftLSpec A B (si : ShiftL A B) '{Equiv A} '{Equiv B} '{RingOne A} 
'{RingPlus A} '{RingMult A} '{RingZero B} '{RingOne B} '{RingPlus B} ;= { 
shiftLproper : Proper ((=) (=) => (=)) (<) ; 
shiftLO :> Rightldentity «) ; 
shiftLS :Vxn,x<(1+n) = 2*x<n}. 

We do not specify shifti as shifti x n = x * 2 " n since on the dyadics we cannot take a negative 
power while we can shift by a negative integer. Since theory on shifting with exponents in 
N and Z is similar we want to avoid duplication of theorems and proofs. To this end we 
introduce a class describing the bi-induction principle. 

Class Biinduction R '{Equiv R} '{RingZero R} '{RingOne R} '{RingPlus R} : Prop 
:= biinduction (P: R Prop) '{IProper ((=) iff) P} ; P ^ (V n, P n o P (1 + n)) ^ V n, P n. 

Since this class is inhabited by any integer and natural implementation we can parametrize 
theory on shifti as follows. 

Context '{SemiRing A} '{ILeftCancellation (.*.) (2:A)} '{SemiRing B} '{IBiinduction B} '{IShiftLSpec A B si}. 
Lemma shiftl_base_plus xyn:(x + y)<n = x<n + y<n. 
Global Instance shiftlJnj: V n, Injective (<n). 

4.4. Decision procedures. The Decision type class collects types with a decidable equal- 
ity |SvdWll| . 

Class Decision P := decide: sumbool P (^ P). 

Using this type class we can declare a parameter '{V x y. Decision (x < y)} to describe a decider 
for < and say decide (x < y) to decide whether x < y or not. This type class allows us to easily 
define additional deciders, like the one for the strict order. We have to be careful however. 
Consider the order on the dyadics. 

Global Instance dy_precedes: Le Dyadic := A (x y : Dyadic), 
ZtoO (mant x) * 2 " (expo x) < ZtoQ (mant y) * 2 " (expo y) 

Now, decide (x < y) is actually (©decide Dyadic (x < y) dyadic_dec, where dyadic_dec is the compu- 
tational conclusion of the decision. Due to eager evaluation, and the absence of dead code 
removal, the second argument, x < y, is also evaluated. Evaluation of this argument results 
in a conversion of x and y into Q, as described above. But since this argument is just a 
proposition it is later thrown away. We avoid this problem introducing a A-abstraction. 

Definition decide_rel '(R : relation A) {dec : V x y, Decision (R x y)} 
(x y : A) : Decision (R x y) := dec x y. 

We can now define: 

Context '{FullPseudoOrder A} '{ITrivialApart A} '{V x y, Decision (x < y)}. 
Global Program Instance lt_dec: V x y, Decision (x < y) | 9 := A x y, 

match decide_rel (<) y x with 

I left E right _ 

I right E ^ left _ 

end. 

Notice that this problem also occurs if we had defined the following bool-based variant. 
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Class BoolDecider (P : Prop) := booLdecide: bool. 

Class BoolDecision (P : Prop) '{BoolDecider P} := booLdeclsion_prf : booLdecide P = true o P. 

4.5. Explicit type casts. The Cast type class collects (explicit) type casts. 

Class Cast A B := cast: A ^ B. 
Implicit Arguments cast [[Cast]]. 
Notation "' x" := (cast _ _ x) (at level 20). 
Instance: Params (@cast) 3. 

This definition allows us to refer to a cast from x : A to B by using an apostrophe, or writing 
cast A B X. An example of an instance of this class is: 
Instance NonNegJnject: Cast (R+) R := @proj1_sig R _. 

Here, is the non-negative cone of an ordered ring R. Contrary to COQ's built-in coercion 
mechanism, our type casts are explicit instead of implicit and type classes are used to register 
them. Our approach has some advantages: 

(1) By using type classes to register casts, we are allowed to parametrize classes with 
casts. An example is the AppRationals class, as defined in Section [5] 

(2) Implicit coercions often introduce ambiguity. Since our approach allows us to refer 
to casts by a (canonical) name, e.g. cast B C (cast A B x), we can avoid this ambiguity. 

(3) Casts can be put in partially applied position, e.g. order_preserving (cast Z Q). 
Coq's coercion mechanism does not allow us to define a coercion from R^ to R nor a 

coercion from a ring to its polynomial ring. More generally, it does not allow most forms 
of parametrized coercions nor non-uniform coercions. If COQ would allow parametrized 
coercions like NonNegJnject, it would have to avoid an infinite loop: to type check x : R, we 
try to type check x : R+, but to do so we try x : (R"'')"'', . . . We suffer from such loops if we 
compose our Cast classes automatically as well. Hence we refrain from adding: 
Instance cast_comp_base '{f : Cast A B} : ComposedCast A B := f. 

Instance cast_comp_step '{f : Cast B C} '{g : ComposedCast A B} : ComposedCast A C := A x, f (g x). 

Matita |ASCTZ07 ] allows parametrized coercions and avoids the loop by not applying coer- 
cions recursively, but instead building a well-chosen set of set of composite coercions [ TasQ Sj . 
Non-uniform coercions |STllj are available in Matita. They are implemented using unifi- 
cation hints, a feature similar to type classes. 

5. The real numbers 

To make our implementation of the reals independent of the underlying dense set, 
we provide an abstract specification of approximate rationals inspired by the notion of 
approximate fields which is used in the Haskell implementation of the exact reals by 
Bauer and Kavler [BK09]. In particular, we provide an implementation of this interface by 
dyadics based on COQ's machine integers. 

Our interface describes an ordered ring containing Z that is dense in Q. Here Z are the 
binary integers from Coq's standard library, and Q are the rationals based on these binary 
integers. We do not parametrize by arbitrary integer and rational implementations because 
they are hardly used for computation. 
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Also, for efficient computation, this interface contains the operations: approximate 
division, normahzation, an embedding of Z, absolute value, power by N, shift by Z, and 
decision procedures for both equality and order. 

Class AppDiv AQ := app_div : AQ ^ AQ ^ Z ^ AQ. 
Class AppApprox AQ := app_approx : AQ Z AQ. 

Class AppRationals AQ {e plus mult zero one inv} '{Apart AQ} '{Le AQ} '{Lt AQ} 

{AQtoQ : Cast AQ Q_as_MetricSpace} '{lAppinverse AQtoQ} {ZtoAQ : Cast Z AQ} 
'{lAppDIv AQ} '{lAppApprox AQ} '{lAbs AQ} '{IPow AQ N} '{IShiftL AQ Z} 
'{V X y : AQ, Decision (x = y)} '{V x y : AQ, Decision (x < y)} : Prop := { 

aq_ring :> @Ring AQ e plus mult zero one inv ; 

aqJriviaLapart :> TrivialApart AQ ; 

aq_order_embed :> QrderEmbedding AQtoQ ; 

aq_strict_order_embed :> StrictOrderEmbedding AQtoQ ; 

aq_ring_morphism :> SemiRing.Morphism AQtoQ ; 

aq_dense_embedding :> DenseEmbedding AQtoQ ; 

aq_div : V x y k, ball (2 " k) ('app_div x y k) ('x / 'y) ; 

aq_compress : V x k, ball (2 " k) ('app_approx x k) ('x) ; 

aq_shift ;> ShiftLSpec AQ Z «) ; 

aq_nat_pow :> NatPowSpec AQ N (") ; 

aqJnts.mor :> SemiRing.Morphism ZtoAQ }. 

Following O'Connor | ;O'C07j . we define the real numbers as the completion of the approx- 
imate rationals. To create functions on the real numbers, we use the monadic operations 
bind or map. This approach is convinient because equality and inequality are decidable on 
the approximate rationals, whereas it is not on the real numbers. For binary functions, e.g. 
addition and multiplication, we use the map2 function, as described in [Q'C07j . 

O'Connor |O'C07 ] keeps the size of the rational numbers small to avoid efficiency prob- 
lems. He introduces a function approx x e that yields the 'simplest' rational number between 
x — e and X + e. We modify the approx function slightly: app_approx x k yields an arbitrary 
element between x — 2'° and x + 2*^. Using this function we define the compress operation on 
the real numbers: compress := bind (A x e, app.approx x (Qdlog2 e)) such that compress x = x. 

In Section 5.4 we will explain our choice of using a power of 2 to specify the precision 
of app.div and app_approx. 



5.1. Order and apartness. Following |BB851 IO'C09| . we define non-negativity and the 
order on the real numbers as follows. 

NonNeg x := Ve : (Q+, —e<xe 

X < y := NonNeg {y — x) 

Bishop and Bridges f BB85j define positivity as the dual of non- negativity: 3e : Q+, e < x e. 
O'Connor [O'C09] defines positivity and the strict order differently so as to avoid a poten- 
tially expensive computation, namely x e — e, to obtain a witness between and x. 

Pos X := {e : Q+ \ £ < x} 

X <T y ■= Pos (y — x) 

We use the T subscript to emphasize that the relation lives in Type. Next, we define 
X Xt y ■= X <T 2/ V y <T X. Extraction of a witness e E (0, x] from Pos x allows us to 
define the reciprocal function of type Vx : R, Xt x — )• R. 
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In order to use our type class based hierarchy we need a strict order and apartness 
relation in Prop. We need this restriction because COQ's present implementation of setoid 
rewriting does not allow rewriting in Type-based relations (see Section 4.1). Our definition 



is similar to Bishop and Bridges', but defined in such a way that we can easily prove a 
correspondence with O'Connor's. 

a; < y := 3n : N, 1 < -n < (y - (1 < {-n - 1)) 

X X y := X < y y y < X 



Using constructive indefinite description (see Section 2.2), it is an easy job to prove that 
we indeed have x < y -H- x <t y and x X y -fr^ x Xx y- Similar to O'Connor |O'C09| . we 
implement a tactic that automatically proves strict inequalities. The tactic terminates iff 
the inequality holds and has quite some similarities with our use of linear search to obtain 
X <T y from x < y. 



5.2. Implementation using the dyadics. The dyadic rationals are numbers of the shape 
n*2'^ for n, e G Z. In order to remain independent of an integers implementation, we abstract 
over it. For our eventual implementation of the approximate rationals we use COQ's machine 
integers, bigZ. Now given an arbitrary integer implementation Int it is straightforward to 
define the dyadics. Here we will just show the ring operations. 
Notation "x \ p" := (exist _ x p) (at level 20). 
Record Dyadic := dyadic { mant : Int ; expo : Int }. 
Infix "$" := dyadic (at level 80). 

Global Instance dyJnject: Cast Int Dyadic := A x, x $ 0. 
Global Instance dy_opp: Grouplnv Dyadic := A x, -mant x $ expo x. 
Global Instance dy_mult: RingMult Dyadic := A x y, mant x * mant y $ expo x + expo y. 
Global Instance dy_0: RingZero Dyadic := cast Int Dyadic 0. 
Global Instance dy_1 : RingOne Dyadic := cast Int Dyadic 1 . 
Global Program Instance dy_plus: RingPlus Dyadic := A x y, 
if decide_rel (<) (expo x) (expo y) 

then mant x + mant y < (expo y - expo x) t _ $ min (expo x) (expo y) 
else mant x < (expo x - expo y) t _ + mant y $ min (expo x) (expo y). 

In this code shifti has type Int lnt+— ^ Int, where Int"'" is a S-type describing the non-negative 
cone of Int. Therefore, in the definition of dy_plus we have to equip expo y — expo x with a 
proof that it is in fact non-negative. 



5.3. Implementation using the rationals. Our development contains additional imple- 
mentations of the AppRationals class using Coq's old rational numbers Q and the new rational 
numbers bigQ (which are built from the machine integers bigZ) . Although creating these im- 
plementations is uninteresting from a performance point of view, it confirms that it is trivial 
to change the underlying dense set from which our real numbers are built. 

To implement the app_approx function in an efficient manner, we use shifts on the under- 
lying integers. Furthermore, to keep the size of the results of the division operation small, 
we incorporate the app.approx function. 
Instance bigQ.div: AppDiv bigQ := A x y, app.approx (x / y). 
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5.4. Power series. Elementary transcendental functions as exp, sin, In and atan can be 
defined by their power series. If the coefficients of a power series are alternating, decreasing 
and have limit 0, then we obtain a fast converging sequence with an easy termination proof. 
For -1 < X < 0, 

Ex 
J 

i=0 

is of this form. To approximate exp x with error e we take the partial sum until ^ < e. 
In order to implement this efficiently we use a stream representing the series and define a 
function that sums the required number of elements. For example, the series 1, a, a^, a^, ... 
is defined by the following stream. 

CoFixpoint powers_help (c : A) : Stream A := Cons c (powers.help (c * a)). 
Definition powers : Stream A := powers.lielp 1 . 

Streams in COQ, like lists in Haskell, are lazy. So, in the example the multiplications are 
accumulated. 

Since COQ only allows structural recursion (and guarded co-recursion) it requires some 
work to convince COQ that our algorithm terminates. Intuitively, one would describe the 
limit as an upperbound of the required number of elements using the Exists predicate. 

Inductive Exists A (P : Stream A Prop) (x : Stream) : Prop := 
I Here : P x — ^ Exists P x 
I Further : Exists P (tl x) — ^ Exists P x. 

This approach leads to performance problems. The upperbound, encoded in unary format, 
may become very large while generally only a few terms are necessary. Due to vm.compute's 
eager evaluation scheme, this unary number will be computed before summing the series. 
Instead O'Connor [O'C09] uses LazyExists. 

Inductive LazyExists A (P : Stream A Prop) (x : Stream A) : Prop := 
I LazyHere : P x — > LazyExists P x 

I LazyFurther : (unit LazyExists P (tl x)) LazyExists P x. 

Unfortunately, our experiments showed that the above still yields too much overhead due 
unnecessary to reduction of proofs. To remedy this issue we introduce the following function 
where Str_nth_tl n s takes the n-th tail of the stream s. 

Fixpoint LazyExistsJnc '{P : Stream A Prop} 

(n : nat) s : LazyExists P (Str_nth_tl n s) LazyExists P s := 
match n return LazyExists P (Str_nth_tl n s) — >• LazyExists P s with 

I O ^ A X, X 

I S n => A ex, LazyFurther (A LazyExistsJnc n (tl s) ex) 
end. 

This function adds n additional LazyFurther constructors. When instantiated with a big 
enough n, computation will suffer from the implementation limits of COQ (e.g. a stack 
overflow) or runs out of memory, before it ever refers to the actual proof. Using LazyExists.inc 
we are able to compute on average twice the amount of decimals as we did before on examples 
such as the ones in Table [2j 

O'Connor's InfiniteAlternatingSum s returns the real number represented by the infinite 
alternating sum over s, where the stream s is decreasing, non- negative and has limit 0. 
We extend this in two ways. First, we generalize various notions to abstract structures. 
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Secondly, as we do not have exact division on approximate rationals, we extend the algo- 
rithm to work with approximate division. The latter requires changing InfiniteAlternatingSum s 
to InfiniteAlternatingSum n d which computes the infinite alternating sum of the stream Ai, ^. 
This allows us to postpone divisions. Also, we have to determine both the length of the 
partial sum and the required precision of the divisions. To do so we find a k such that: 

B| (app_div Hk dk (log^) + ^) 0. (5.1) 

Now k is the length of the partial sum, and ^ is the required precision of division. Using 
O'Connor's results we have verified that these values are correct and such a k indeed exists 
for a decreasing, non- negative stream with limit 0. 

As noted in Section [5} we have specified the precision of division in powers of 2 instead 
of using a rational value. This allows us to replace (5.1) with: 

B| (app_div Hk dk (log e - {k + 1)) + 1 <^ (log e - {k + 1))) 0. 

Here k is the length of the partial sum, and 2\ where / = log e — (A; + 1), is the required 
precision of division. This variant can be implemented without any arithmetic on the 
rationals and is thus much more efficient. 

This method gives us a fast way to compute the infinite alternating sum, in practice, 
only a few extra terms have to be computed and due to the approximate division the 
auxiliary results are kept as small as possible. 

Similarly, using this method to compute infinite alternating sums, we use the following 
series to implement atan x and sin x for x G [—1, 1]. 

atan x = y~^(— 1)' 



1=0 

oo 



X 

s\n x = > (— 1)' * 



(2i + l)! 

2i+l 



2i + l 



i=0 

We extend these functions to their complete domain by repeatedly applying the following 
formulas [O'CJCQj . 

exp X = (exp (x < 1))^ (5.2) 

exp X = r (5-3) 

exp (— x) 

X ( x\3 

sin X = 3 * sin — — 4 * (^sin —J (5-4) 

atan x = —atan (— x) (5-5) 
vr 1 

atan x = atan — for < x (5-6) 

2 X 

vr /X — 1\ 

atan x = — — atan for < x (5.7 

4 Vx + 1/ ^ ' 

Since we do not have exact division on the approximate rationals, we parameterize infinite 



sums by two streams in Equation 5.4, 5.6 and 5.7 



The serie s described in this section converge faster for arguments closer to 0. We use 
Equation 5.2 and 5.4 repeatedly to reduce the input to a value |x| G [0,2^^). For 50 < fc. 



this yields nearly always major performance improvements, for higher precisions setting it 
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to 75 < yields even better results. Unfortunately, we are unaware of a similar trick for 
atan. We define tt in terms of atan using the following Machin-like formula. 

Ill 1 
TT := 176 * atan— + 28 * atan—— — 48 * atan—— + 96 * atan: 



57 239 682 12943 

Again, here we notice the purpose of parameterizing infinite sums by two streams. We 
define cos in terms of sin. 

cos X = 1 — 2 * ^sin —j 

O'Connor |O'C07t rO'C09j subtracts multiples of 27r to reduce the arguments of sin and cos. 
In our tests this did not lead to performance improvements because our implementation of 
vr turned out to be slower than the performed range reductions. 

5.5. Square root. We use Wolfram's algorithm [Wol021 p. 913] for computing the square 
root. Its complexity is linear, in fact it provides a new binary digit in each step. 

Context '(Pa : 1 < a < 4). 
Fixpoint AQrootJoop (n : nat) : AQ * AQ := 
matcli n with 
I O ^ (a, 0) 
I Sn ^ 

let (r, s) := AQroot_loop n in 
if decide_rel (<) (s + 1) r 
then ((r - (s + 1 )) « {2:Z), (s + 2)< (1 :Z)) 
else (r « (2:Z), s « (1 :Z)) 
end. 

Let us write (r„, Sn) for the n-th pair of approximations. We prove the following facts by 
induction: 

si + 4r„ = r+^a (5.8) 
r„ < 2s„ + 4 (5.9) 
2"^s„<s„+„<2™(s„ + 4)-4 (5.10) 
r„ < 2^+" (5.11) 



By |5.8| (2 (" +^)sn )^ + 2 ^"r„ = a. By |5.1l| 2 ^"r„ converges to as n tends to oo. 



Therefore, by 5.10, 2^("+^)s„ is a Cauchy sequence which furthermore converges to the 
root of a. 

Next, we extend the square root to its entire domain by repeatedly applying the fol- 
lowing formula. 

I~x 



Vx = 2* ^, ^ 
"4 

O'Connor's COQ implementation |O'C08j includes the much faster Newton iteration, 
whose complexity is logarithmic in the number of decimals. The function to iterate is: 
Definition f (r : Q) : Q := r / 2 + a / (2 * r). 

Because of the absence of exact division on our approximate rationals we cannot implement 
this function directly. However, we can implement it on our real numbers. As the above 
definition does not use sharing, we have to define this function on the reals by first defining: 
Definition f (r : AQ) (e : Qpos) : AQ := (r + approx.div (Qdlog2 e) a r) < (-1). 
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POl 


^'\n ( <^in ( ^'\n 1 i i 


5.000 


25s 


2.3s 


P02 




5.000 


3.3s 


1.7s 


P03 


sin e 


5.000 


13s 


1.2s 


P04 


exp (vr * \/l63) 


5.000 


22s 


2.0s 


P05 


exp (exp e) 


5.000 


43s 


2.6s 


P06 


log (1 + log (1 + log (1 + log (1 + tt)))) 


500 


107s 


2.5s 


P07 


exp 1000 


20.000 


1.1s 


0.7s 


P08 


cos (10^0) 


20.000 


6.7s 


1.4s 


P09 


. /o log 640320 ^ 


5.000 


33s 


16s 


Pll 


tan e + atan e + tanh e + atanh ^ 


500 


41s 


3.2s 


P12 


asin - + cosh e + asinh e 


500 


99s 


3.2s 



Table 1. Haskell, compiled with ghc version 6.12.1, using -02. 





Expression 


Decimals 


Old 


New 


Decimals 


New 


POl 


sin (sin (sin 1)) 


25 


46s 


0.6s 


500 


3.8s 


P02 


\/vr 


25 


0.3s 


0.03s 


500 


6.8s 


P03 


sin e 


25 


36s 


0.1s 


500 


1.9s 


P04 


exp (vr * \/l63) 


10 


214s 


0.1s 


500 


3.7s 


P05 


exp (exp e) 


10 


36s 


0.2s 


500 


3.2s 


P07 


exp 1000 


10 


2662s 


1.0s 


2.000 


4.9s 


P08 


cos (10^0) 


25 


lis 


0.3s 


2.000 


12s 



Table 2. COQ trunk, revision 14023. 



and then showing that it gives rise to a continuous function / : AQ — t- AR which we finally 
lift to a function bind / : AR — t- AR on the reals. In this way we take care of sharing, division 
and intermediate use of the approx function (see Section [s]) all in one go. We hope the future 
correctness proof to be quite smooth, since we work with exact real numbers. We have 
implemented this in HASKELL, which performs really well. 

6. Benchmarks 

The first step in this research was to create a Haskell prototype based on O'Connor's 
implementation of the real numbers in HASKELL |O'C07j . The second step was to implement 
and verify this prototype in COQ. Our COQ development contains verified versions of: the 
field operations, exponentiation by a natural, computation of power series, exp, atan, sin, 
cos, vr and the square root. 

In this section we present some benchmarks, taken from the 'Many Digits' friendly 
competition |NW09j . comparing the old and the new implementation, both in Haskell 
and COQ. All benchmarks have been carried out on an Intel Core Quad 2.4 GHz with 
8GB of memory running Debian GNU/LiNUX with kernel 2.6.38. The sources of our 
developments can be found at http : //robbertkrebbers . nl/ research/reals. 
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Table [T] shows some benchmarks in Haskell with compiler optimizations enabled (-02) 
and Table [2] compares our COQ implementation with O'Connor's. More extensive bench- 
marking shows that our HASKELL implementation generally benefits from a 15 times speed 
up while the speed up in COQ is generally more than a 100 times for small examples already. 
This difference between the comparison of the Haskell and the COQ implementation is ex- 
plained by the fact that O'Connor's Haskell implementation already uses rational numbers 
built from fast integers and incorporates various optimizations, while his COQ implemen- 
tation does not. The last column of Table [2] indicates that our new implementation is able 
to compute an order of magnitude more decimals in the same amount of time. 

We also compared the new reals built from Coq's fast rationals (Section 5.3) and our 
dyadic rationals (Section 5.2). For exp, sin and cos we obtain quite similar results due to the 
our range reductions to reduce the length of the power series. In case of the square root, 
the dyadics rationals are much faster because wolfram iteration is designed for an efficient 
shift. It is interesting to notice that tt and atan benefit the least from our improvements, as 
we are unaware of range reductions to reduce the length of the series. 

We conclude this section with a comparison between the performance of Wolfram's 
algorithm in COQ and Haskell. The Haskell prototype (without compiler optimizations) 
is quite fast, computing 10,000 iterations (giving 3,010 decimals) of \f2 takes 0.2s. In COQ 
it takes 7.4s using type classes and 7.2s without type classes. Here we exclude the time 
spend on type class resolution. Thus type classes cause only a 3% performance penalty on 
computations. 

Unfortunately, the COQ implementation is slow compared to Haskell. Laurent Thery 
suggested that this is due to the representation of the fast integers, which uses a tree 
with a fixed depth and when the size of the integer becomes too big uses a less optimal 
representation. Increasing the size of the tree representation and avoiding an inefficiency in 
the implementation of shifts reduces this time to 5.4s. 



7. Conclusions and Related work 

We have greatly improved the performance of real number computation in COQ using 
Coq's new machine integers. We produced highly structured and abstract code using 
type classes with no apparent performance penalty. Moreover, COQ's notation mechanism 
combined with Unicode characters gives nicely readable statements and proofs. Type classes 
were a great help in our work. However, the current implementation of instance resolution 
is still experimental and at times too slow (at compile time). 

Canonical structures provide an alternative, and partially complementary, implemen- 
tation of type classes [GZNDlT] . By choice, canonical structures restrict to deterministic 
proof search, this makes them more efficient, but also somewhat more intricate to use. The 
use of canonical structures by the Ssreflect team [GGMR0"9] makes it plausible that with 
some effort we could have used canonical structures for our work instead. However, the 
SSREFLECT-library is currently not suited for setoids which are crucial to us. The integra- 
tion of unification hints [ARC TO 9] into COQ may allow a tighter integration of type classes 
and canonical structures. 

We needed to adapt our correctness proofs to prevent the virtual machine from eagerly 
evaluating them. Lazy evaluation for Prop would have allowed us to use the original proofs. 
Moreover, setoid rewriting over relations in Type would have made our work much easier. 
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The experimental native_compute performs evaluation by compilation to native OCAML 
code. This approach uses the Ocaml compiler available and is interesting for heavy com- 
pilation. Our first experiments indicate a 10 times speed up with Wolfram iteration. Un- 
fortunately, nativexompute does not work with COQ trunk yet, so we were unable to test it 
with our implementation of the reals. 

The Flocq project [iBMllj formalizes infinitary floating-points in COQ. It provides 
a library of theorems on multi-radix multi-precision arithmetic and supports efficient nu- 
merical computations inside COQ. However, the current library is still too limited for our 
purposes, but in the future it should be possible to show that they form an instance of our 
approximate rationals. This may allow us to gain some speed by taking advantage of fine 
grained algorithms instead of our more straightforward ones. 

The encoding of real numbers as streams of 'bits' is potentially interesting. However, 
currently there is a big difference in performance. The computation of 37 decimals of the 
square root of 1/2 by Newton iteration |JP09j . using the framework described in |Ber07| 
IJulOSj . took 12s. This should be compared with our use of the Wolfram iteration, which 
gives only linear convergence, but with which we nevertheless obtain 3,000 decimals in in a 
similar time. On the other hand, the efficiency of vr in their framework is comparable with 
ours. Berger |Ber09j . too, uses co-induction for exact real computation. 

The present work is part of a larger program to use constructive mathematics based 
on type theory as a programming language for exact analysis. This should culminate in 
a numerical ODE-solver. As part of this goal, we plan to build a type class interface for 
metric spaces. 

Cohen and Mahboubi |CMllj provide a formalization of quantifier elimination for the 
theory of decidable real closed fields. Together with the formalization of cylindrical algebraic 
decomposition |Mah07j this should lead to an efficient decision procedure. Such a procedure 
would have been useful in our current formalization effort. 

Acknowledgements. We thank Eelis van der Weegen for many discussions and Pierre 
Letouzey and Matthieu Sozeau for closing some of our bug reports. 
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